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Instability of dilute granular flows on rough slope 
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We study numerically the stability of granular flow on a rough slope in collisional flow regime in 
the two-dimension. We examine the density dependence of the flowing behavior in low density 
region, and demonstrate that the particle collisions stabilize the flow above a certain density 
in the parameter region where a single particle shows an accelerated behavior. Within this 
parameter regime, however, the uniform flow is only metastable and is shown to be unstable 
against clustering when the particle density is not high enough. 
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Granular flow on a slope is one of the simplest situ- 
ations to see the characteristic behavior of granular dy- 
namics. When the inclination angle is smaller than a 
certain value (the angle of repose), the material never 
flows because of the stress sustained by the friction. Be- 
yond that angle, the surface layers of the material may 
flow like a fluid, but the bottom part of the materials 
may remain solidified when the inclination angle is not 
steep enoughff If the inclination is increased further, all 
of the material starts to flow rapidly, and the interaction 
between particles or between a particle and the slope is 
dominated by the inelastic collision, rather than friction. 

Many researches have been done in such a collisional 
flow regime, but mgs±_Qf them focus on the property 
of the uniform flow.ErtJQ ) In such researches, the depth 
dependence of the flow properties such as velocity or 
density profile is investigated, assuming that the flow 
is uniform in the direction along the slope. It is known, 
however, that the granular materials have the tendency 
to cluster due to the inelastic collision, which causes the 
formation of density waves in the case of granular flow in 
a vertical pipeBEP Therefore, it is natural to expect that 
this tendency will cause some instability in the uniform 
flow on a slope. 

Granular flow of an independent particle, or a single 
particle behavior, has been studied, and has been found 
to show three types of motion dependiru^pn the inclina- 
tion angle and roughness of the slope.Lr&Lp For the fixed 
roughness of the slope, the following behaviors are ob- 
served upon increasing the inclination angle: (i)The par- 
ticle stops after a few collisions with the slope for any ini- 
tial velocity (regime A). (ii)The particle quickly reaches 
a constant averaged velocity in the direction along the 
slope and shows almost steady motion; the averaged ve- 
locity does not depend on the initial condition (regime 
B). (iii)The particle jumps and accelerates as it goes 
down the slope (regime C). 

In the present work, we study how the above single 
particle picture is modified in the collisional flow with 
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finite density by the particle collisions. Based on nu- 
merical simulations, we determine the parameter region 
where the uniform collisional flow is realized with the 
finite density, and examine the stability of the uniform 
flow. _ ._. 

We employed the discrete element methodliliElP with 
the normal and the tangential elastic force and dissipa- 
tion. In the simulations, granular particles are modeled 
by two-dimensional disks with mass m and diameter d. 
When the two disks i and j at positions r j and rj with 
velocities Vi and Vj and angular velocities Ui and uij are 
in contact, the force acting on the particle i from the 
particle j is calculated as follows: The normal velocity 
v n and the tangential velocity Vt are given by 

v n = v i3 • n, v t = Vij ■ t - (d/2)(uJi + uij), (1) 

with the normal vector n — rjj/|r,j| = (n x ,n y ) and 
the tangential vector t = (—n y ,n x ). Here, r t j = Tj — 
r, and Vij = Vj — Uj. Then the normal force and 
the tangential force Ffj acting on the particle i from the 
particle j are given by 

F% = -k n (d - \r tj \) + r) n v n , (2) 

4 =mm(|/i t |,/x|i^|)sign(/ lt ), (3) 

with 

h t = k t u t + r] t v t , Ut = j v t dt. (4) 

Here, Ut is the tangential displacement, and to is the time 
when the particles started to contact. k n and k t are elas- 
tic constants, r\ n and r\t are damping parameters, and \i is 
the Coulomb friction coefficient for sliding friction. Each 
particle is also subject to the gravity, and the gravita- 
tional acceleration is given by g = g(sin8, — cos 9). The 
surface of the slope is made rough by gluing the parti- 
cles identical with the moving ones with spacing 2ed with 
e = 0.001 (Fig. 1). The periodic boundary condition is 
imposed in the x direction. The system size L is deter- 
mined by the number of the particles glued on the slope 
n s as L = (1 + 2e)dn s . 

All quantities which appear in the following are given 



2 



Namiko Mitarai and Hiizu Nakanishi 




Fig. 1. A schematic illustration of the system geometry. 



Table I. Parameters used in the simulations. 

k n k t T] n r\ t p 

5 x 10 4 (5 x 10 4 )/3 35.68 35.68/3 0.5 



in the non-dimensionalized form in terms of the length 
unit d, the mass unit to, and the time unit r = (d/g) 1 / 2 . 
The parameters used in the simulations are tabulated in 
Table I. With these parameters, the normal restitution 
coefficient (the ratio of normal relative velocities of pre- 
and post-collision) e„ turns out to be around 0.7. We 
integrated the equations of motions using the predictor- 
corrector method in the third order with a constant time 
step 5t = lx 1CT 4 . 

We first investigate how the uniform collisional flow 
appears when we increase the density of flowing particles 
from the one particle limit. 

As we have referred, a particle on a slope shows three 
types of behavior depending on the value of 9. It is 
known, however, that the boundary between the regime 
B (the constant velocity regime) and C (the acceleration 
regime) is not sharp; there is a parameter region where 
the particle can attain steady state or can accelerate, 
depending on the initial condition. We should also note 
that, when the initial kinetic energy of the particle is too 
small, the particle may stop in the regime B and in part 
of the regime C. In the present model with the param- 
eters in Table I, it turns out that the regime A roughly 
corresponds to the region sin 9 < 0.11, the regime B to 
0.11 < sin0 < 0.14, and the regime C to 0.16 < sinfl. 

Now we focus on what happens when the number of 
moving particles n is increased with the finite L. It is 
easy to expect that the dissipation becomes larger upon 
increasing n, because the collisions between particles or a 
particle and the slope become more frequent. Therefore 
it is obvious that all the particles will stop irrespective 
of the initial condition if 9 is in the regime A. Even for 
any 9 within the regime B, it turns out that the collisions 
between particles eventually stop all the particles. 

We examine closely the flowing behavior for the larger 
inclination angle 9 in the regime C. The situation is pre- 
sented by showing the simulations with sin# = 0.45 and 
L = 20.04. We set the initial condition as follows: 

Xi(0) = (L/n)(i - 1), y< (0) - (1 + 2e)d + o&, (5) 
Ui(0) = 0, Vi (0) = 0, wj(0) = 0, (6) 
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Fig. 2. The semi-log plot of E(t) for (A)p = 0.50, (B)p = 0.60, 
(C)p = 0.65, and (D)p = 1.00, with sin6» = 0.45 and L = 20.04. 
E(t) continues to grow when p < 0.60, while it is bounded when 
p > 0.65. 



where Xi(t) and Uj(i) («»(*)) are the coordinate 

and the velocity of center of mass of the ith particle at 
time t in the x (y) direction, respectively. is a ran- 
dom number uniformly distributed over the interval [0, 
1]. The appropriate value of a is a few times the diam- 
eter, and we adopted a = 3 ~ 4 in actual simulations. 
In the following, we define the 'density' of the particles 
as p = n/L. In the simulations with L = 20.04, we do 
not observe non-uniformity along the slope in the parti- 
cle distribution, therefore the parameter p is enough to 
describe the situation. In order to characterize the qual- 
itative difference of accelerated behavior in the regime C 
and the uniform flow, we investigate the p dependence 
of the averaged kinetic energy E(t) defined as 

n 

£(*) = -£Si(*), (7) 

where Ei{t) is the kinetic energy of the ith particle at 
time t; 

Ei{t) = l -m [ui{tf + v t {t) 2 } + ^Jwi(t) 2 , (8) 

with I = md 2 /8. 

In Fig. 2, the time evolutions of E(t) are shown for 
several values of p. When p is small (p < 0.60), E(t) 
grows rapidly and continuously, but the growth rate be- 
come smaller as p increases. In this region, each par- 
ticle jumps and rarely collides with each other. When 
p > 0.65, E(t) still grows rapidly in early stage, but 
its long-time behavior seems to be bounded and fluctu- 
ate around a constant value. In this case each particle 
also jumps, but often collides with other particles and is 
prevented from jumping up infinitely. Based on this ob- 
servation, we define the uniform flow in the low-density 
limit as the flow which is uniform along the flow direction 
with the value of E(t) being bounded, namely the energy 
dissipation due to inelastic collisions balances with the 
energy gain from the gravity. 

The average value of the kinetic energy E in the uni- 
form flow should be related to X, the distance between 
collisions along the slope, by (1 — efyE ~ mgXsin9, be- 
cause the energy gain by the gravitation should balance 
with the loss due to inelastic collisions. The simulation 
shows the typical value of X to be 0(1) ~ O(10) in the 
uniform flow, thus we expect E to be 0(1) ~ O(10). On 
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Fig. 3. The 'phase diagram' obtained from the simulations with 
L = 20.04. The accelerated regime (open box), the uniform flow 
regime (filled box), and the regime where all the particles stop 
(plus) arc shown when more than 50 % of trials result in the 
corresponding behavior. 



the other hand, E diverges in the accelerated regime. 

Figure 3 shows the 'phase diagram' summarizing the 
behavior in terms of 9 and p obtained from numerical 
simulations. In order to obtain the diagram, 10 simula- 
tions with L = 20.04 were done for each 9 and p. Wc 
identify the three behaviors, namely (i) accelerated mo- 
tion, (ii) uniform flow, and (iii) static state where all 
the particles come to rest, by the following criterion; the 
system is determined to be in the accelerated motion (in 
the static state) if E(t) exceeds 500 (becomes zero) by 
t = 20, 000. It is determined to be the uniform flow if 
neither of them happen. The final behavior sometimes 
depends on the initial condition. This, we believe, is due 
to the finite size effect. Each point in Fig. 3 is deter- 
mined to be in one of these regimes if more than 50 % of 
the trials show the one of the above behaviors. It is ex- 
pected that, if we could take the L — > oo limit with fixed 
p with keeping the particle distribution uniform along 
the slope, such initial condition dependence should dis- 
appear. 

Now we examine the stability of the uniform flow by 
taking the system size large enough to observe non- 
uniformity along the slope. First we show the simula- 
tion results of the flow in the system with sin# = 0.45, 
L = 1002 and n = 1000, i.e. p = 1.0. The initial con- 
dition is also given by eq. (6), therefore the particles 
are uniformly distributed along the slope at first. Figure 
4 shows the time evolution of the flux <&(i) defined as 
the number of the particles which pass x — L/2 during 
the time interval At — 10. After a short transient time, 
the flux becomes almost constant with small fluctuations 
(300 < t < 2000), which indicates uniform flow is real- 
ized. Then $(t) shows large fluctuation which grows in 
the course of time (2000 < t < 3000); this means the 
clustering behavior is triggered by the fluctuation. Fi- 
nally $(£) begins to oscillate almost periodically with 
large amplitude (t > 3000). This oscillation of the flux 
indicates that one large cluster of particles travels in the 
system with almost constant velocity; we can see that the 
cluster is stable once it is formed. From this observation, 
we expect that the uniform flow is not stable against clus- 
tering. Considering the fact that the uniform flow has a 
finite life time during which the flow seems fairly stable, 
we expect that a fluctuation of finite size is necessary to 



Fig. 4. The time evolution of the flux &(t) in the simulation with 
L = 1002 and p = 1.0. A fluctuation grows to form a cluster 
spontaneously. 
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Fig. 5. The time evolution of (a) E(t) and (b) E\ (t) in the sim- 
ulation with L = 1002 and p = 1.0. Compared with Fig. 4, 
we can see that, after the clustering occurs (t > 3000), the time 
averaged value of E(t) becomes larger than that in the uniform 
flow stage, and Ei(t) begins to show large fluctuation. 



trigger the clustering, which means the uniform flow is 
metastable. 

In Figs. 5, the time evolutions of (a) the averaged 
kinetic energy E(t) and (b) the kinetic energy of a par- 
ticular particle Ei (t) are also shown. Both of them are 
almost constant in the stage of the uniform flow. Then 
E(t) increases in 2000 < t < 3000. In t > 3000, the fluc- 
tuation of Ei(t) becomes considerably large, and E(t) 
begins to fluctuate around another constant value which 
is larger than the one in the stage of the uniform flow. 
The reason why the kinetic energy increases when a clus- 
ter is formed can be understood as follows: When the 
cluster is formed, the region with the density lower than 
the threshold value to prevent the acceleration (<~ 0.65 
in the case of sin# = 0.45, see Fig. 3.) appears locally, 
and particles in such a spatial region can be highly ac- 
celerated. However, the particle will be caught in the 
cluster sooner or later, and then quickly lose its kinetic 
energy. This mechanism maintains the moving cluster, 
and results in the large fluctuations in E\ (t) . 

also examined the system size and the density depen- 
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Fig. 6. The time evolution of the flux <t>(t) in the simulation with 
L = 501 and (a) p = 1.0, (b) p = 2.0, and (c) p = 0.75. 



dence of the clustering. The time evolutions of the flux 
for a few values of densities with L = 501 are also shown 
in Figs. 6. Figure 6 (a) shows the result with p = 1.0 
(n = 500), which is the same density as the simulation in 
Fig. 4. In Fig. 6 (a), however, no clustering behavior can 
be seen clearly even though the fluctuation is very large. 
We have also simulated in the system with L — 250.5 
and p = 1.0 (n — 250), but the clustering behavior was 
not found, either. On the other hand, the density de- 
pendence of the system behavior can be seen in Figs. 6 
(b) (p = 2.0) and (c) (p = 0.75). It is found that the 
uniform flow is maintained in Fig. 6 (b), while the clear 
oscillation of the flux associated with the cluster forma- 
tion is seen in Fig. 6 (c). From these results, we can see 
the general tendency that the uniform flow tends to be 
stabilized as the system size is smaller or as the density 
is higher. 

In summary, we have examined the two-dimensional 
granular flow on a rough slope in the collisional flow 
regime by numerical simulations. It was shown that the 
mutual collisions among particles stabilizes the flow even 
in the accelerated regime for a single particle system. 
The phase diagram was determined for the accelerated, 
uniform flow, and stopping regime in terms of the par- 
ticle density p and the inclination angle 9 for a partic- 
ular system size. The stability of uniform flow was also 
examined and we found that a large single cluster ap- 
pears spontaneously out of a uniform initial state. It 
was shown that the smaller the particle density is, or the 



larger the system size is, the less stable the uniform flow 
is. 

Before concluding, let us make a comment on the sys- 
tem size dependence of the instability. Consider the sit- 
uation where one large cluster, whose length along the 
slope is I and typical hight is h, is formed. If we can 
neglect the particles outside the cluster, we can assume, 
as the first approximation, that the volume of the clus- 
ter (oc Ih for 2-dimensional case) is proportional to n. 
Namely we get the relation (l/d) (h/d) oc n = (pd)(L/d). 
Therefore, if (l/d) and (h/d) are scaled with (L/d) for 
fixed p as (l/d) oc {L/d) 13 and (h/d) oc (L/d) 13 , we get 
the relation f3' = 1 — (3. Because both (l/d) and (h/d) 
should be increasing function of (L /d) for fixed p, we get 
an inequality for the exponent (3 that < f3 < 1. Then 
there should exist an critical system size L c such that 
l(L c ) ~ L c ; a cluster can exist only in the case of L > L c 
for which l(L) < L, and the uniform flow become stable 
when L < L c . Using similar argument, we can also pre- 
dict the existence of the critical density p c for a fixed L. 
However, we should note that the moving cluster seems 
to be maintained by the balance between the outgoing 
flux of the particles at the front of the cluster and incom- 
ing flux at the tail, therefore the motion of the particles 
outside the cluster may affect on the shape of the clus- 
ter in the periodic boundary condition. Therefore, the 
discussion here may not simply hold and the exponent 
13 may depend on p. The scaling behaviors should be 
confirmed more carefully. 

The spontaneous cluster formation out of uniform 
flow is seen also in the granular flow through ^ verti- 
cal pipetltr and in the traffic flow on a freewayOO It 
is an interesting problem to find out whether those phe- 
nomena and the cluster formation in surface flow have 
a common mathematical structure at phenomenological 
level. 

Part of the computation in this work has been done us- 
ing the facilities of the Supercomputer Center, Institute 
for Solid State Physics, University of Tokyo. 
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